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Summary 

A study of plane-acoustic-wave propagation in small tubes 
with a cross section in the shape of a flattened oval is described. 
To begin, theoretical descriptions of a plane wave propagating 
in a tube with circular cross section and between a pair of 
infinite parallel plates, including viscous and thermal damping, 
are expressed in similar form. For a wide range of useful duct 
sizes, the propagation constant (whose real and imaginary parts 
are the amplitude attenuation rate and the wave number, 
respectively) is very nearly the same function of frequency 
for both cases if the radius of the circular tube is the same 
as the distance between the parallel plates. This suggests that 
either a circular-cross-section model or a flat-plate model can 
be used to calculate wave propagation in flat-oval tubing, or 
any other shape tubing, if its size is expressed in terms of an 
equivalent radius, given by g = 2 X (cross-sectional area)/ 
(length of perimeter). Measurements of the frequency response 
of two sections of flat-oval tubing agree with calculations based 
on this idea. Flat-plate formulas are derived, the use of 
transmission-line matrices for calculations of plane waves in 
compound systems of ducts is described, and examples of com- 
puter programs written to carry out the calculations are shown. 


Introduction 

“The problem of the propagation of sound waves in gases 
contained in cylindrical tubes is a classical one, to which 
famous names are connected like Helmholtz, Kirchhoff, and 
Rayleigh. Since then many papers have been written on the 
subject, often in relation to studies dealing with the dynamic 
responses of pressure transmission lines.” So states H. 
Tijdeman, referring to 34 previous papers in the first paragraph 
of one of his reports on the subject (ref. 1). (See also refs. 
2 and 3.) 

Much less has been written about the propagation of sound 
waves in ducts of noncircular cross section, however, 
particularly for propagation in narrow tubes, where viscous 
and thermal losses at the walls are significant. It was the goal 
of this study to develop an analytical framework for calculating 
the pressure-transfer properties of the flat-oval tubing shown 
in cross section in figure 1, so that its behavior could be 
predicted when it is used as probe tubing for a pressure 


transducer. The hope was that results could be cast in a form 
similar to the recursion relation of Bergh and Tijdeman (ref. 
4), which has been used to analyze pressure probes of circular 
cross section (ref. 5). This recursion relation is obtained in 
a slightly more general form in this report (appendix C). 

However, in the section Solving Fluid Equations for 
Noncircular Tubing, it is argued that a Bergh and Tijdeman 
type recursion relation cannot be found for tubing of 
noncircular cross section for a viscous, thermally conducting 
gas. That is to say, the fundamental differential equations 
describing the behavior of a fluid in a duct cannot be solved 
in closed form to the same degree of approximation used for 
the Bergh and Tijdeman analysis in any coordinate system 
other than circular cylindrical. 

The solution for a plane wave propagating between two 
infinite parallel plates can be found, however, and the section 
Comparison With Flat-Plate Geometry is concerned with a 
comparison of the circular-tube and parallel-plate results. Over 
a useful range of tubing sizes and frequencies, the two solutions 
are remarkably similar if a generalized radius g = 2 x (area 
of cross section)/(length of perimeter of cross section) is 
properly used to characterize the size of the waveguides. Thus, 
within the range where the solutions are similar, either should 
give a close approximation to the acoustical properties of flat- 
oval tubing, described by its generalized radius. This con- 
clusion applies also to tubing with other cross-sectional shapes 
as well. 

The pressure-transfer properties of short sections of flat- 
oval tubing were measured, and the results are presented in 
the section Measurements on Flat-Oval Tubing. The theoretical 
predictions agree very well with the measurements. Limits to 
this technique for calculating the pressure-transfer properties 
of noncircular tubing are suggested in the section Concluding 
Remarks. 

Four appendixes conclude the report. A list of symbols used 
in the equations is given as appendix A. Appendix B gives 
details of the solution for plane waves propagating between 
infinite parallel plates. Appendix C describes the use of 
transmission-line matrices for calculating the acoustical 
properties of systems of tubing. (This method is equivalent 
to but more straightforward and powerful than the Bergh and 
Tijdeman recursion relation.) Listings of several computer 
programs written to implement the transmission-line 
calculations are presented and explained in appendix D, with 
several numerical examples for reference. 


l 



Figure 1.— Cross section of flat-oval tubing investigated. 


Solving Fluid Equations for 
Noncircular Tubing 

There are four basic equations describing the properties of 
a viscous gas with a finite thermal conductivity. These four 
equations are the following: 

The force equation, which includes viscosity, is the Navier- 
Stokes equation (refs. 6 and 7), which may be written 


, Du' 


(V • V)u' +-V(V *u') (1) 


In this equation, which is a vector equation having three 
components, p' is the fluid density, u' is the fluid velocity, 
p' is the absolute pressure, and p is the constant shear viscosity 
coefficient. The Stokes assumption that the bulk (or volume) 
viscosity is zero has been made, and the differential operator 


D d 

— = — + u' • V 
Dt dt 


( 2 ) 


In curvilinear coordinates, the scalar operator V • V can be 
evaluated by using the identity 


(V • V)u' = V(V -u') - VX(VXU') (3) 


The equation of continuity of the fluid (or conservation of 
mass) may be written (ref. 8, p. 8, p. 512) 


Dp' 

— — h p' V • u' =0 
Dt 


(4) 


The thermodynamic properties of the gas are assumed to 
be described by the ideal gas equation of state 


p' =p'RoT' 


(5) 


and an energy equation describing heat conduction (refs. 9 and 
10 ) 


P c p 


DT ' Dp' 
Dt Dt 


- x v 2 r =o 


( 6 ) 


Here T' is absolute temperature, c p is the specific heat at 
constant pressure, and X is the thermal conductivity of the 
fluid; a second-order term describing heat transfer due to 
internal friction has been neglected (see refs. 1 and 4). 

When written in cylindrical coordinates, equations (1), (4), 
(5), and (6) form the starting point for the analysis in the 
appendixes of references 1 and 4, where it is assumed from 
symmetry that there is no azimuthal velocity. The solutions 
for the circular tube proceed by setting 


p' = Pa + pe jost T' = T a + Te J(at 

p' = p a + pe jo3t u' = ue** 

where j = V- l, and p a , p a , and T a are constants represent- 
ing the ambient or average values of pressure, density, and 
temperature. As usual, the position-dependent disturbance 
amplitudes p, p, and T are assumed to be small compared to 
the ambient values. Similarly, the magnitude of the fluid 
velocity amplitude u is assumed small compared to c, the free- 
space adiabatic phase velocity of sound given by c-\lyp a /p~, 
where y is the ratio of specific heats y = c p /c v . In addition, 
it is assumed that (1) the internal tube radius r is small 
compared to the free-space wavelength (a >rlc « 1), (2) the 
radial velocity component is smaller than the axial velocity 
by about this same factor, and (3) the flow is laminar 
throughout the system. (Laminar flow is interpreted to mean 
that differentiation with respect to the axial coordinate yields 
a quantity of the same order as does dividing by the free-space 
wavelength; differentiation with respect to the radial coordinate 
gives a result of the same order as does dividing by r.) Under 
these assumptions, closed-form solutions to these equations 
in cylindrical coordinates are obtained in references 1 and 4. 

For purposes of this investigation, it is natural to consider 
solving the basic equations to the same degree of approxi- 
mation in other coordinate systems. The difficulty encountered 
in such attempts may be illustrated by using rectangular 
coordinates as an example. 

Choose the x-axis to lie along the center axis of a long tube 
with rectangular cross section, with the y- and z - axes 
perpendicular to the duct walls. Denote the x-component of 
the harmonic amplitude of the fluid velocity by u, the y- 
component by v, and the z-component by w. Writing the basic 
equations (1), (4), (5), and (6) in rectangular coordinates, 
substituting equations (7), and retaining only the most 
significant terms in each equation (as in the circular tube case) 
yield the following six equations: 
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1 dp p ( d 2 u d 2 u\ 

J0}U JL + JL + 

Pa fa Pa\dy fa / 

(8) 

O 

II 

1 

^ 1^3 

(9a) 

O 

II 

(9b) 

/ du dv dw\ 
jup - -p a y— + j y + y z j 

(10) 

p = -^(p - p<A>t) 

(11) 

( d 2 T d 2 T\ 

j*^~ x W + &) +Jap 

(12) 


These are similar to equations (6) through (10) in the 
appendix of reference 4 or equations (Bl) through (B5) of 
reference 1, except that they contain two independent co- 
ordinates or components perpendicular to the tube axis instead 
of one. 

To this approximation as in the circular case, equations (9a) 
and (9b), which are the two transverse components of the 
Navier-Stokes equation, directly imply that the pressure waves 
are plane waves. The planes of constant pressure coincide with 
planes of constant x. However, the rectangular case differs 
from the circular case because it has less symmetry. 

The rectangular case has two independent transverse 
components of fluid velocity, denoted by v and w, whereas 
by symmetry the circular case has only one transverse velocity, 
the radial velocity. In order to solve the problem, then, we 
need one more equation in the rectangular case than in the 
circular geometry to determine an additional unknown. This 
additional equation is the third component of the Navier-Stokes 
equation given in equation (9), but unfortunately it provides 
no information about v or w. The two transverse velocities 
appear only together in equation (10) with no additional 
relation which can be used to separate them. In addition, if 
one tries to solve equation (12) for temperature fluctuations, 
one finds that with both y- and z-derivatives present the 
separation of variables procedure is not successful. 

Thus, the equations do not appear solvable to this 
approximation in rectangular coordinates. Furthermore, we 
would expect the same difficulty in any other coordinate system 
lacking circular symmetry and thus having two independent 
transverse velocity components, if we seek a plane wave in 
pressure. In spite of the Bessel functions which appear, the 
circular geometry is fundamentally the simplest case to 
consider and is the only shape duct for which a closed solution 
is possible. 


Comparison with Flat-Plate Geometry 

As mentioned by Lord Rayleigh, another effectively two- 
dimensional case which can be solved, in addition to the 
circular duct, is a plane wave propagating between a pair of 
infinite parallel plates (ref. 11). A solution for that case is 
derived in appendix B, and the results may be summarized 
as follows: If a plane pressure wave propagates in the jc- 
direction between a pair of isothermal, rigid, infinite planes 
located at y = ±/i/2, the ^-dependent amplitude of the pressure 
disturbance at frequency is 

p(x) = Ae* x + Be~** (13) 

Here, A and B are complex constants determined by boundary 
conditions, and the complex propagation constant <p is 

(14) 

In this expression, 



and Pr is the Prandtl number. (See appendix B.) At any 
particular x, the average amplitude of the longitudinal fluid 
velocity is 

u(x) = — F(a) Ue** - Be -**) (18) 

juPa v 7 

Interestingly enough these expressions are of precisely the 
same form as the corresponding quantities in the circular duct 
solution, as found in equations (29), (30), and (42) in the 
appendix of Bergh and Tijdeman (ref. 4). The only difference 
is that for the infinite-flat-plate geometry the function F has 
replaced the Bessel function ratio J 2 /Jq appearing in the 
circular-tube case. The arguments of the respective functions 
are the same if one chooses the spacing h between the plates 
equal to the radius r of the circular tube. 

Not only are the overall expressions of the same form, but, 
for equal arguments of the type used in these equations, the 
magnitude and phase of the function F are actually quite similar 
to the corresponding quantities for the Bessel function ratio. 
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If both are evaluated for an argument written as £=j y2 s, 
where s is real (this is the particular form of a), then 


.. m js 2 

(19) 


(20) 

lim — = lim F(f) = - 1 + - 

S ~°%(D s ~°° t 

(21) 


For large arguments, both have the limiting value — 1 , so that 
for very wide conduits, very high frequencies, high ambient 
pressure, or vanishing viscosity and thermal conductivity the 
propagation constant becomes yco/c . Between the extreme 
limits, the magnitudes and phases of the two quantities are 
plotted in figures 2 and 3. There it can be seen that the two 
quantities are not identical, but are surprisingly similar. To 
indicate the frequency range in which the important variation 
of these two functions occurs in the system of interest, figures 
2 and 3 s how th e values of s appropriate for 100 and 329 Hz 
if s = hyJp a o>/fi for air at atmospheric pressure and 300 °C 
and h = 0.061 cm. 

Quantities that are calculated by using these two similar 
functions are also very much alike. Graphs of the real and 
imaginary parts of the propagation constant <p as a function 
of frequency are shown in figures 4 through 6, for several 



Figure 2.— Comparison of magnitudes of J 2 (^)/J 0 (^) 311(1 [(tan f/2)/(£/2)] — 1 
as function of real parameter s , where f = j m s. 



Figure 3.— Comparison of phases of J 2 ( £)/Jo ( and [(tan $V2)/(f/2)] - 1 as 

function of real parameter s , where f — j m s. 

values of r = h 9 given in centimeters. The imaginary part of 
<p is the wave number, inversely proportional to the phase 
velocity of the wave in the conduit. It is the primary quantity 
determining the resonant frequencies of standing waves in a 
tube. The real part of <p is the factor giving the amplitude 
attenuation rate with distance. In figures 4 through 6, the solid 
straight line is the wave number for an ideal, inviscid gas, 
the dashed lines ending in a circle are computed by using the 
circular-tube expressions, and the dashed lines ending in a 
square show the infinite-flat-plate values. The most noticeable 
difference between the circular-tube and flat-plate propagation 
constants is that for the smallest tubing the flat-plate geometry 
has a larger attenuation factor, or greater viscous and thermal 
damping, than the circular case. But the smallest case shown 
has such large damping that tubing of that size would not be 
useful for pressure-transmission purposes. The intermediate 
pair of graphs turns out to show the propagation constant 
appropriate for the flat-oval tubing of interest, and there the 
flat-plate damping is only slightly greater than the circular- 
tube damping. 

For the particular case of infinite flat plates, figures 4 
through 6 illustrate the assertion that theoretical expressions 
derived for waves propagating in a circular tube may be 
adapted to describe the properties of waves in ducts of arbitrary 
cross-sectional shape. This may be done by simply replacing 
r, the radius of the circular tube, by a generalized quantity 
calculated as 

g — 2 X (cross-sectional area of duct)/(length of 

perimeter of cross section) (22) 

wherever r appears in the equations (refs. 12 through 14). For 
a circular tube, g = r, while in the limit that the width of a 
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(a) Wave number for circular-cross-section tube with radius of 1 .0 cm, wave 
between flat plates separated by 1.0 cm, and plane wave without thermal 
or viscous damping. 

(b) Attenuation factor for circular-cross-section tube with radius of 1.0 cm, 
and wave between flat plates separated by 1 .0 cm. 

Figure 4.— Wave number and attenuation factor, imaginary and real parts, 
respectively, of propagation constant of equations (13) and (14) as functions 
of frequency. 


rectangular duct becomes infinite, g = h. The assertion seems 
true at least in that p(x) and u(x) are very nearly the same 
function of g , frequency, and gas properties for waves 
propagating both in a circular duct and between infinite parallel 
plates over a wide range of sizes. It would be reasonable to 
expect p(x ) and u{ x) to depend in the same way on the same 
parameters in flat-oval tubing or tubing of other cross-sectional 
shapes as well. 



Frequency, Hz 

(a) Wave number for circular-cross-section tube with radius of 0.060 cm, 
wave between flat plates separated by 0.060 cm, and plane wave without 
thermal or viscous damping. 

(b) Attenuation factor for circular-cross-section tube with radius of 0.060 cm, 
and wave between flat plates separated by 0.060 cm. 

Figure 5.— Wave number and attenuation factor, imaginary and real parts, 
respectively, of propagation constant of equations (13) and (14) as functions 
of frequency. 


Shown in figure 7 are rough sketches of circular, flat-oval, 
and infinite-plate waveguides drawn with the same value of 
g. The pressure and average-velocity waves in these three 
systems, then, are described by almost the same expressions. 
The properties of systems made up of tubing with any cross- 
sectional shape should be calculable with a recursion relation 
like that of Bergh and Tijdeman (ref. 4) or equivalent 
transmission-line formulas (see appendix C), where in most 
places r is replaced by the appropriate g. 

However, care must be exercised, for not all such 
replacements would be correct. The exceptions occur in 
expressions for the characteristic impedance of a tube found 
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Frequency, Hz 

(a) Wave number for circular-cross-section tube with radius of 0.010 cm, 
wave between flat plates separated by 0.010 cm, and plane wave without 
thermal or viscous damping. 

(b) Attenuation factor for circular-cross-section tube with radius of 0.010 cm, 
and wave between flat plates separated by 0.010 cm. 

Figure 6.— Wave number and attenuation factor, imaginary and real parts, 
respectively, of propagation constant of equations (13) and (14) as functions 
of frequency. 

from equations (CIO), (C12), and (C19) or the V t tube 
volume factors in the Bergh and Tijdeman recursion relation 
given here as equation (C57), where essentially the cross- 
sectional area S of the tube enters as a result of requiring 
conservation of mass flowing through the tube. The mass per 
unit time passing point x is given by 


— = P a Su(x) (23) 

at 


In adapting this requirement to tubing of different shape, even 
if 5 is expressed as xr 2 , set S equal to the actual cross- 
sectional area of the new tube, rather than changing r to g. 
Similarly, when calculating an appropriate end correction to 
the geometrical length of a tube (0.82 r for flanged tubing or 



Figure 7.— Cross-sectional shapes of circular, flat-oval, and infinite-parallel- 
plate ducts, all having same size generalized or equivalent radius. 


0.61 r for an unflanged opening), replace r by Vs7r for a 
noncircular tube (ref. 8, pp. 348-350). 

Measurements on Flat-Oval Tubing 

In order to evaluate the applicability of the circular-tube 
and/or flat-plate expressions to the case of flat-oval tubing, 
two samples of such tubing were mounted on the resonator 
chamber of the large Galton whistle sound source shown in 
figure 8 and described in reference 15. Silicon strain-gauge 
pressure transducers were used to measure the pressure in the 
whistle and in a small cavity at the end of the tubing. The 
geometrical length of one sample was 5.08 cm and of the other 



Figure 8.— Galton whistle sound source for testing frequency response of flat- 
oval probe microphone tubing. 
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was 7.62 cm, while the cavity volumes were 0.01325 and 
0.01630 cm 3 , respectively. The cross-sectional area of the 
tubing was 0.01714 cm 2 , and the length of the perimeter was 
0.571 cm. Thus, g — 0.060. When the expected pressure- 
transfer properties of the tubing were calculated, an end 
correction of 0.07 cm was added to the geometrical length of 
each tube. 

The magnitudes of the pressure disturbance amplitudes were 
measured by using two Princeton Applied Research model 126 
lockin amplifiers operated as ac voltmeters, and the ratio of 
the end-cavity pressure amplitude to the whistle pressure 
amplitude was plotted as a function of frequency. The relative 
phase of the two signals was also measured with one of the 
lockin amplifiers. Care was taken to subtract the phase shift 
in the signal channel input filter of the amplifier at each 
frequency. 

The experimental data are shown as crosses superimposed 
on theoretical predictions in figures 9 (shorter tube) and 10 
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(a) Ratio of amplitude of transducer cavity pressure to amplitude of whistle 
pressure at probe tube inlet. 

(b) Relative phase of transducer cavity pressure and whistle pressure at probe 
tube inlet. 


Figure 9.— Ratio of cavity to whistle pressure amplitude and relative phase 
of two signals as functions of frequency for flat-oval tubing of geometrical 
length 5.08 cm and equivalent radius g = 0.060. 
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(a) Ratio of amplitude of transducer cavity pressure to amplitude of whistle 
pressure at probe tube inlet. 

(b) Relative phase of transducer cavity pressure and whistle pressure at probe 
tube inlet. 

Figure 10.— Ratio of cavity to whistle pressure amplitude and relative phase 
of two signals as functions of frequency for flat-oval tubing of geometrical 
length 7.62 cm and equivalent radius g - 0.060. 


(longer tube). The lengths of the bars on the crosses in these 
figures do not signify experimental uncertainty. There are, in 
principle, two theoretical curves on each graph, one from the 
circular tube formulas and one for infinite flat plates, both 
evaluated for g = 0.060 cm. The two curves are essentially 
the same, however, except near the peaks in the pressure- 
amplitude ratios. There the slightly greater damping in the flat- 
plate case leads to a slightly lower peak. The relative phase 
curves are totally indistinguishable. 

The agreement between the experimental data and the 
theoretical curves is excellent. Since figure 7 suggests that the 
flat-oval tubing may resemble a pair of flat plates more than 
a circular tube, one might expect the experimental points to 
lie closer to the curves of the flat-plate model than those of 
the circular one. In fact they do, but the possibility of slight 
extra, uncontrolled damping or air leakage in the apparatus 
implies that such fine distinctions would be hard to defend. 


Concluding Remarks 

Very satisfactory predictions of the acoustical properties of 
pressure-transmission systems incorporating flat-oval tubing 
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Figure 11.— Calculated magnitude of ratio of transducer cavity pressure to 
pressure at probe tube inlet as a function of frequency for flat-oval tubing 
of length 7.7 cm and equivalent radius g — 0.020 cm. 


with an equivalent radius g = 0.060 cm can be made by using 
equations derived either for circular tubes or a pair of infinite 
flat plates. This result was anticipated from the surprising 
similarity of the theoretically calculated propagation constants 
for waves in circular tubes and waves between flat plates, if 
both types of ducts possess the same equivalent radius as 
defined in equation (22). Experimental measurements of the 
pressure-transfer properties of samples of flat-oval tubing with 
g = 0.060 over a frequency range from 500 to 3500 Hz agreed 
very well with the theoretical calculations, confirming the 
applicability of the theoretical models investigated. 

Further calculations show, however, that for only somewhat 
smaller values of g the greater damping in the flat-plate model 
leads to significant differences between the flat-plate and 
circular-tube predictions. Figures 11 and 12 display predicted 
pressure-amplitude ratios for 7.7-cm lengths of tubing with 
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Figure 12.— Calculated magnitude of ratio of transducer cavity pressure to 
pressure at probe tube inlet as a function of frequency for flat-oval tubing 
of length 7.7 cm and equivalent radius g = 0.010 cm. 


g = 0.020 and 0.010 cm, respectively. The damping is higher, 
and predictions from the two models are significantly different. 
This suggests a lower bound of about g = 0.040 cm, below 
which circular-tube properties may no longer be adapted with 
confidence to tubing of other cross-sectional shapes with the 
equivalent radius formalism. Further experiments would be 
required to determine if either of the models considered here 
gives a good representation of very small flat-oval tubing. 


National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio, March 13, 1986 
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Appendix A 
Symbols 


A,B complex constants of integration (eq. (13)) 

C U C 2 complex constants of integration (eq. (B13)) 

C capacitance per unit length of transmission line 
(eq. (Cl)) 

c adiabatic phase velocity of sound (eq. (11)) 

c p specific heat at constant pressure (eq. (6)) 

c v specific heat at constant volume 

D/Dt differential operator (eq. (2)) 

F geometry-dependent function for flat plates (eq. 

(16)) or ratio of Bessel functions for circular tube 
(eq. (CIO)) 

F^x) arbitrary function of integration (eq. (B27)) 
f(x) function for separation of variables (eq. (BIO)) 

8 shunt leakage conductance per unit length of 
transmission line (eq. (C2) 

g equivalent radius (generalized radius) (eq. (22)) 
h distance separating parallel plates (eq. (17)) 

i position-dependent harmonic amplitude of V (eq. 

(C4)) 

i' electric current in transmission line (eq. (C6)) 
Jq 9 J 2 integer-order Bessel function of first kind (eq. 
(19)) 

j imaginary unit 

K complex constant of integration (eq. (Bll)) 

k polytropic constant for cavity (eq. (C39)) 

L length of section of tube or duct (eq. (C23)) 

£ inductance per unit length of transmission line 
(eq. (C6)) 

f subscript labeling tube section (eq. (C46)) 

m mass of parcel of gas (eq. (23)) 

n complex effective polytropic factor for tube (eq. 

(B31)) 

Pr Prandtl number (eq. (B6)) 

p position-dependent harmonic amplitude of excess 

pressure (eq. (7)) 

p a average absolute pressure (eq. (7)) 

p f absolute pressure (eq. (1)) 

q(z ) function for separation of variables (eq. (BIO)) 

q r electric charge (eq. (Cl)) 

R electric resistance 

R 0 ideal gas constant (eq. (5)) 


(R electric resistance per unit length of transmission 
line (eq. (C6)) 

r inside radius of tube with circular cross section 

S cross-sectional area of tube (eq. (23)) 
s general real parameter (eq. (19)) 

T position-dependent harmonic amplitude of 

temperature variations (eq. (7)) 

T a average or ambient absolute temperature (eq. (7)) 
T absolute temperature (eq. (5)) 

t time (eq. (7)) 

U harmonic volume velocity amplitude (flow 

amplitude) (eq. (C8)) 

u axial component of u (eq. (8)) 

u average value of u over cross section of tube (eq. 
(18)) 

u position-dependent harmonic amplitude of fluid 

velocity (eq. (7)) 
u' fluid velocity (eq. (1)) 

V position-dependent harmonic amplitude of V' (eq. 
(C4)) 

V volume of cavity (eq. (C39)) 

V t volume of tube section (eq. (C56)) 

V' electric potential difference (eq. (Cl)) 

v,w transverse rectangular components of u (eq. (10)) 

x axial rectangular coordinate (eq. (8)) 

total shunt electric or acoustic admittance per unit 
length (eqs. (C5) and (02)) 

y,z transverse rectangular coordinates (eq. (8)) 

Z(x) acoustic impedance of tube at x (eq. (Cl 7)) 

Z c characteristic impedance of tube (eq. (09)) 

Z in acoustic input impedance of tube (eq. (C32)) 

Z T acoustic terminating impedance of tube (eq. 

(C33)) 

% total series electric or acoustic impedance per 
unit length (eqs. (C7) and (CIO)) 

z dimensionless scaled transverse coordinate (eqs. 

(B8) and (B20)) 

a shear wave number (eq. (B7)) 

7 ratio of specific heats, c p lc v (eq. (11)) 
f general complex argument of function (eq. (16)) 

A thermal conductivity (eq. (6)) 
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Appendix B 

Plane Waves Propagating Between Infinite Planes 


Consider a rectangular coordinate system with two infinite 
rigid plates located at y = ±hl2. Suppose a wave propagates 
in the space between them, and choose the jc-axis to point in 
the direction of propagation. The linearized equations de- 
scribing the gas are those given as equations (8) through (12) 
in the section Solving Fluid Equations for Noncircular Tubing, 
except that since any ^-location is equivalent to any other, the 
z-derivative terms and hence the z-component fluid velocity 
amplitude w disappear from the equations. Thus we seek to 
find the small deviations from equilbrium by solving 


This may be solved by separation of variables. Let 

T(x, z) =f(x)q(z) (BIO) 

Substituting this into equation (B9) yields 


d 2 q 1 p{x) 


<7(Z)+Ti = — 


dZ 2 P a C p fix) 


= K 


(Bll) 


with the conclusions 


1 dp p d 2 u 

JUU= ~7 a to + 7ad? 

0=-^ 

dy 

(Bl) 


fix) = „Pix) 

Pa c pK 

(B12) 

(B2) 

and 



( du dv\ 

JUf>= - p \y x + 7y) 

(B3) 


q(z) — Ci sin z + C 2 cos z + K 

(B13) 

II 

i 

h 

(B4) 

If the plates are isothermal, so that the temperature fluctations 
vanish at the walls, the two boundary' conditions are 

„ , & 2 T • 

jo>p a c p T= X—^+jup 

(B5) 


T= 0 

(B14) 


We start with equation (B5). When the Prandtl number 
which relates viscous and thermal effects is defined as 


for 


Pr.& 


<B6) 


or 


the shear wave number as 


for 




PaU 


(B7) 



<7 = 0 


z = 


orfPr 

± 

2 


(B15) 


and a new dimensionless variable z as 


Thus we find that in equation (B13) 


z = 



equation (B5) may be written as 


(B8) 



(B16) 


„ d 2 T p 
T+ — r = — 
dz P a Cp 


yielding 

(B9) 


ll 


T(x,z) =f(x)q(z) 


Notice that at any particular x, the average longitudinal fluid 
velocity is 


1 - 


cos z 



pM 

PaCp 


(B17) 



'hi 2 


u(x,y)dy 


In terms of the original variables, this is 


T(x t y) = 


1 - 



P(x) 

Pa c p 


1 

WPa 



(B23) 


(B18) 

The quantity in square brackets appears in several expressions, 
and it is given a special function notation: 


Inserting the result of equation (B18) in equation (B4) and 
using the fact that c p — c v = Rq , we find 



(B24) 


p(*,y) = 




Finally we look at the equation of continuity, equation (B3), 
which may be written 


X 


1 



dv p du 

dy ^ p a dx 


(B25) 


(B19) Inserting equation (B19) for p into equation (B25) and equation 
(B22) for u, we get 


Next, the Navier-Stokes equation, equation (Bl), may be 
solved by separation of variables as was the energy equation. 
If we let 


z = 



equation (Bl) becomes 

d 2 u 1 dp 

u + — z - — ~ 

dz JUPa “X 


(B20) 


(B21) 



Letting u(x,z) =f(x)q(z) and requiring that u vanish at the 
plates because of viscous drag, we find 


u(x,y) = : 

JUPa 


cos 


( a y\ 

\T) 


r a ' 
COS I - 


- 1 


dp 

dx 


(B22) 



If we then integrate with respect to y 9 the result is 
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p(x) = Ae <fiX + Be~ ,fiX 


(B29) 


v(*,30 =— ( y- 

JUPa\ C 


1 - 1 

y 



where the propagation constant <p may be written 



by using the notation 


(B30) 



At the rigid plates, v, the transverse component of the fluid 
velocity, must vanish. Setting equation (B27) equal to zero 
first at y = +/i/2 and then at y = -h!2 and adding the two 
results show that F^jc) = 0. Subtracting one equation from 
the other and cancelling h give 


F(a) TT ~—y- 

dx z c 




p = 0 (B28) 



(B31) 


The quantity n may be regarded as a complex effective 
poly tropic factor for pressure changes in the tube (ref. 4). 

Since <p is a complex constant with positive real part, 
equation (B29) together with equation (7) represents a 
superposition of damped traveling plane waves, one of 
amplitude A traveling in the -jc-direction and one of amplitude 
B traveling in the +*-direction. The constants A and B are 
determined by boundary conditions. 

Inserting equation (B29) and its derivative into the 
expressions for v (eq. (B27», u (eq. (B22)), p (Eq. (B19», 
and T (eq. (B18)) completes the solution of the problem. 


This has the solution 
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Appendix C 

Transmission-Line Formalism 


Once the general expressions for the pressure wave in a tube 
have been obtained, it remains to calculate the transfer 
properties of particular arrangements of probe tubing and 
transducer cavities. Bergh and Tijdeman (ref. 4) derived a 
workable but cumbersome recursion formula for this purpose. 
Another approach, which is fairly commonly used for one- 
dimensional acoustic waves propagating in ducts, is based on 
equations developed to describe electrical transmission lines. 
Since the transmission-line formalism with its ABCD matrices 
is much more compact and flexible than the Bergh and 
Tijdeman recursion relation of reference 4 (the matrices enable 
one to easily calculate pressure transfer relations for a greater 
variety of combinations of tubing and to evaluate other 
quantities of interest such as input impedances), it seems 
worthwhile to show how the acoustical equations are analogous 
to the transmission-line ones. The acoustical analysis may then 
be carried forward in a form suggested by the formalism 
developed for transmission lines. 


Equations for Electrical Transmission Lines 

Consider a coaxial cable transmission line for definiteness, 
although the same concepts apply to parallel-wire and other 
types of transmission line (ref. 16). Any transmission line is 
characterized by a distributed electric capacitance per unit 
length ( 3 , inductance per unit length £, resistance per unit 
length (R (including both conductors), and leakage conductance 
per unit length g of the gas or other dielectric between the 
conductors (g Ax is the reciprocal of the inter-conductor 
leakage-resistance in length Ax). If a potential difference 
between the two conductors is set up at one end of the cable, 
charge will start to flow. Let us find differential equations for 
the current i ' in the central wire and potential V' of the central 
wire relative to the outer one, as functions of the distance x 
along the wire and of the time. 

As suggested in reference 16, choose an element of the 
transmission line of infinitesimal length Ax for analysis. There 
are four physical relations or quantities of interest: the 
capacitive relation between the voltage V' and the charge on 
the central wire in Ax, the inductive emf associated with the 
rate of change of the current in Ax, the i'R drop in potential 
along the conductors, and the leakage current between the 
conductors due to conductance through the intervening 
medium. 

The charge q' on length Ax of the central conductor is 

q' = (VAx) V' (Cl) 

The rate at which charge leaves Ax, traveling down the 
conductor, differs from the rate at which it enters by an amount 


equal to the leakage current through the dielectric plus the rate 
at which the local capacitance is being charged. That is, 

Ai' = — V'(QAx) - ~~ (C2) 

ot 

Substituting from equation (Cl) and dividing by Ax give 

ar dv' 

— = -SK' -e— (C3) 

dx at 

If harmonic time dependence is assumed in the form 

V ' = Ve*"* (C4a) 

V = i^ x (C4b) 

then the first basic transmission-line relation between complex 
amplitudes is 

-yv (C5) 

dx 

Here y is the total distributed shunt electric admittance per 
unit length of line. 

The potential drop along Ax is 

ar 

AV' = - i ' ((R Ax) - (£Ax) — (C6) 

Dividing through by Ax and assuming harmonic time de- 
pendence result in the other transmission-line relation: 

dV 

— = — (<R + jo)£)i m - Zi (Cl) 

dx 

In this equation, Z is the total distributed series electric 
impedance per unit length of the line. 


Analogous Quantities for Acoustic Tubes 

In the natural analogy with an acoustical duct, the excess 
pressure p is associated with the electric potential and the 
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volume velocity is associated with the electric current. At any 
cross section the volume velocity or flow V is defined as 

U = Su (C8) 

where S is the cross-sectional area of the duct carrying the 
wave, and u is the fluid velocity averaged over the cross 
section. Multiplying the volume velocity by the average fluid 
density gives the mass passing any point per unit time, so that 
conservation of mass can be conveniently expressed directly 
in terms of conservation of volume velocity. 

Notice from equations (B23) and (B24) that 


dp = jupg v 
dx SF(cc) 


(C9) 


so by analogy with equation (C7), the total distributed series 
acoustic impedance per unit length of the duct may be said to be 


Z = 


. AA PaC _J_ 
J \C ) S F(a) 


(CIO) 


As a function of a, Z depends on the viscosity but not on the 
thermal conductivity of the gas. The imaginary part of Z/co 
is the distributed inertance of the fluid per unit length, while 
the real part of Z is the distributed series acoustic resistance 
of the duct per unit length, due to viscous damping. In the 
section Comparison with Flat-Plate Geometry, it is pointed 
out that the expressions for pressure and longitudinal fluid 
velocity waves in a circular duct deduced by Bergh and 
Tijdeman in reference 4 have the same form as the ones derived 
here for parallel flat plates. All the circular duct results may 
be obtained from the ones for parallel flat plates by replacing 
the function F(a) defined in equation (16) by the Bessel 
function ratio J 2 (a)/J 0 (a). For that reason, it is possible to 
use the equations of this appendix to analyze either flat plates 
or circular ducts, by letting F(ct) refer to the function of 
equation (16) or to the Bessel function ratio, depending on 
the case of interest. In addition, experiments described in this 
report show that, within limits, tubes with other cross-sectional 
shapes may be treated with either flat-plate or circular-tube 
functions if the tube is described in terms of its equivalent 
radius. 

Differenting equation (C9) with respect to x and using 
equation (B28), we find 


dU 

dx 



7 - 1 


FlayfPr 


\P (Cll) 


which is the other transmission-line relation for the acoustical 
wave in a duct. The total distributed-shunt acoustic admittance 
per unit length is 



which depends on thermal conductivity but not viscosity. The 
imaginary part of Tj/w is the distributed compliance of the fluid 
per unit length, while the real part of ^ is the reciprocal of 
a resistance describing energy loss to the walls through thermal 
conduction per unit length. 

The two quantities Z and y define the fundamental physical 
properties of the duct. When multiplied by the appropriate 
length Ax, they give the series impedance and shunt admittance 
of an infinitesimal element of tube -a tube which is very short 
compared to the wavelength of any acoustical disturbance in 
the tube. Since an elastic fluid in a long tube sustains wave 
motion in which a harmonic disturbance at one point is not 
in phase with the disturbance at other locations, Z and y 
cannot be used to directly find the impedances of a finite length 
of tube. Rather, the two quantities Z and ^ are first used to 
calculate the characteristic impedance and propagation constant 
of the tube. It is then possible to analyze finite tubes with 
various terminations, or combinations of tubes. 

We may proceed as follows: Starting from the transmission- 
line type equations 


dp 

dx 


= -ZU 


dU 

dx 


= -yp 


(Cl 3a) 
(Cl 3b) 


we differentiate equation (Cl 3a) and substitute equation (Cl 3b) 
to get 


j3=-z^ = (z<y ) P (Ci4) 

dx 1 dx 

This wave equation for p implies the result of equation (B29), 
with propagation constant 

v = sfzy (ci 5) 

From equations (C13a) and (C13b) again, the corresponding 
volume velocity wave is 


—(Ac** - Be~* x ) 
Zdx Z' f 

=V1(— ) 


(Cl 6) 
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The combination of Z and ^ appearing in the last expression 
also has physical significance. The acoustic impedance of the 
duct at any point is defined as 


Po = A + B 


(B-A) 


z ( ,).£^ 

m 


Given our expressions for the waves. 


-VI 


Be-** + Ae* x 


V <\j Be~ <px — Ae* x 

If a wave is introduced atx - 0 into a semi-infinite tube along 
the positive x-axis there will be no reflected wave traveling 
back toward the origin. This implies that A = 0. The 
impedance of the duct under these conditions is independent 
of position and is called the characteristic impedance of the 
duct Z r . Thus 


p L = Ae* L + Be-* 


(Be~* L - Ae* L ) 


Equations (C23) and (C24) may be solved for A and B as 


(Pl ~ Z c U L )e~* L 


(Pl + Z c U L )e* L 


-Ji-z-z 

y v * y 


Inserting equations (C25) and (C26) into equations (C21) and 
(C22) yields 


At any x , the relation between Z(x) and the characteristic 
impedance of the duct Z c at that point indicates the relative 
amplitude and phase of leftward- and rightward-traveling 
pressure waves in the duct. Rearranging equation (Cl 8) shows 
specifically that the ratio of leftward to rightward wave 
amplitudes is 


Po = Pl cosh <pL + U L Z C sinh <pL (C27) 


U 0 = p L ~ sinh <pL -I- U L cosh <pL (C28) 


which may be written as a matrix equation 


Ae**eP* Z(x) - Z c 
Be~ ipx ^ t ~ Z(x) +Z C 


cosh ipL Z c sinh <pL p L 


Transmission-Line Matrices 

The stipulation of the values of p and its derivative dp /cbc 
at some point determines a unique solution to the wave 
equation, as the constants A and B in equation (B29) are thereby 
determined. Since the volume velocity U is proportional to 
dp/dx (see eq. (Cl 3)), the choice of pressure and volume 
velocity at any point in a duct determines the pressure and 
volume velocity at any other point. That is, the pressure and 
volume velocity at any point in a duct can always be expressed 
as a function which is a linear combination of the pressure 
and volume velocity at a reference point. This fact leads to 
the transmission-line equations and ABCD matrices (ref. 17). 

For example, let the wave amplitudes at x = 0 in a long duct 
be called p 0 and U 0 , with the corresponding quantities at 
x — L being p L and U L . From equations (B29) and (Cl 6) we 
obtain 


U 0 — sinh <pL cosh <pL U L 


The matrix is called the transmission-line matrix, or the ABCD 
matrix, the latter name coming from expressing the generic 
2 by 2 matrix as 


There is no connection between the name ABCD matrix and 
the A and B in equations (C25) and (C26). 

The transmission-line matrix can easily be used to calculate 
the input impedance of a duct of length L which is closed at 
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the far end. The boundary condition is expressed by setting 
U L = 0 . Then 


and expressing the flow into the other branches in a similar 
way, we obtain 


P 0 


U 0 



cosh <pL Z c sinh <pL 


1 


— sinh <pL 


L Z 


cosh <pL 


Pl 


U u = U d +(±r + ±r + ^ + ...)p d (C37) 

vZi Z2 z 3 


or 


Pi cosh (pL 


Pu " 


1 

o' 


- Pd ’ 

— sinh <pL 

l Z c J 

(C31) 

Uu 


1 1 1 

1 1 f ... 

Z\ z> z. 

1 


U d 


(C38) 


and 


Zi„ = 77 = Z c coth <pL 
U 0 


(C32) 


The input impedance can be calculated just as easily for an 
arbitrary terminating impedance Z r . The result is 


Z in = Z ( 


Z T cosh <pL + Z c sinh <pL 
Z T sinh <pL 4- Z c cosh <pL 


(C33) 


Other events in the life of a transmission line may also be 
expressed in terms of ABCD matrices including the effects of 
various kind of discontinuties in the line. Suppose that at some 
point in a duct there are several side branches so that not all 
the mass arriving at that point from the upstream duct leaves 
by the downstream duct. Using subscripts u to describe the 
situation just upstream of the discontinuity, d for the 
downstream side of the discontinuity, and 1 , 2 , 3 , ... for the 
entrances to the side branches, we have 


As pointed out again in the section Applications of ABCD 
Matrices to Compound Systems, this matrix is equal to the 
product of several matrices of the same form, one for each 
of the individual side branches. 

A particular side branch could be an actual side tube attached 
to the main duct to monitor the average total pressure, with 
input impedance calculated as above or with other A BCD 
matrices; or the side branch could be a local cavity, some of 
the flow into which pumps up the pressure rather than going 
out the other side. Bergh and Tijdeman (ref. 4) show, in the 
developments leading up to their appendix equation (53), for 
a cavity of volume V containing gas with polytropic constant 
k (value between 1 for an isothermal process and y for an 
adiabatic one) and ddd (dimensionless diaphragm deflection) 
factor a that the relation between the rate at which mass is 
flowing in and out of the volume and the excess pressure is 


dm jo)y r ( V 

— z=z J —~ V\ o + - }p 

dt c 2 V kj 


(C39) 


Since dm/dt can be considered equal to p a U , the effective 
input impedance of such a cavity is 


and 


Pu - Pd - P\ ~ P2 = Pi = 


(34) 


Z = — = 
z,n U 


A c\yV( 1 

j ~ — (" + 7 
\c/p a c\ k 


(C40) 


U u — Uj + U\ 4- U 2 + t / 3 + ... 


(C35) an d the appropriate ABCD matrix for such a volume is 


Expressing U x in terms of the input impedance to branch 1, 
denoted Z u as 


7 _ Pl _ P d 


(C36) 


/o>\yV/ 1 

/( - — or + - 
\c/p a c\ 


(C41) 
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Many other ABCD matrices can be derived. One can imagine 
a discontinuity in the form of a porous plug in a duct or an 
impervious, flexible membrane with mass, where there is a 
pressure discontinuity but the volume flow is continuous. This 
amounts to a lumped series impedance Z' in the duct, with 
matrix given by 


'ir 

o 1 


(C42) 


There is actually a small series impedance associated with a 
step change in the radius of a circular duct (ref. 18), which 
can be placed in such a matrix. 

An arbitrary discontinuity can be represented by a com- 
bination of series and shunt impedances. Ah ABCD matrix has 
been written for a section of duct with a conically tapered cross 
section (ref. 17). 


Applications of ABCD Matrices to Compound Systems 

Suppose one wishes to analyze the pressure-transfer 
properties or calculate the input impedance of a system of many 
interconnected parts, such as the general one shown in 
figure 13. All that is required is to write down the ABCD 
matrix for each separate part and then multiply them together. 
Since 



Figure 13.— General compound linear acoustic system. 
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a compound side branch discontinuity such as a cavity with 
an attached side tube can be represented by either one overall 
matrix or the product of two matrices, one for each of the 
branches at that point. The pressure and volume velocity at 
the left entrance to tube 2 in figure 13 are given by 


" Pi ' 


matrix for " 


” matrix for " 


' A ' 

- Ui_ 


tube 2 


cavity 1 


0 


(C44) 


while at the far left entrance to the complete system we have 


Po 

U 0 


tube 

8 


hole 

7 


tube 

6 


tube 

5 


side 

tube 

4 


cavity 

3 


ABCD matrices 


It is quite easy to write computer programs to evaluate the 
matrices, multiply them, and print out input impedance 
Zj n = Po/Uo or pressure-transfer function p x /p 0 . (Since the 
equations are linear, for these calculations p x may be set 
equal to an arbitrary constant such as 1.) 


pm d 


tube 

2 


cavity 

1 



P i 
0 


(C 45) 
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Figure 14.— Acoustic system consisting of cavity with input and output tubes, 
used in derivation of Bergh and Tijdeman recursion relation. 


Derivation of Bergh and Tijdeman Recursion Relation 

While it is much easier to program and use the ABCD 
matrices than the Bergh and Tijdeman recursion relation (ref. 
4), it is a satisfying exercise to derive the recursion relation 
from the matrix formalism and thus show they are equivalent. 
Consider the system of two tubes and a cavity shown in 
figure 14. At the entrance to the tube of length L f+1 , just 
after the cavity of volume V b we have from equation (C29) 


Pt 


cosh \Lf+ j 

Zc,t + 1 si* 1 <Pf+i£f+i 

U, 


sinh <pi+ x Lt+ x 

Z c ,t + 1 

cosh <pg +x L (+x 


X 


Pf+i 

U (+ 1 


(C46) 
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Multiplying this equation out, we find from the first component 
Pi = Pe + 1 cos h Vi+iLe+x + ZcjhiUm sinh ^f+i^t+i (C47) 
which implies 


U» 1 = 


Pi ~ Pi + 1 cos h V^+A+i 
Z cf+ i sinh <p M L e+l 


(C48) 


and from the second component 


Uf = sinh <pt+\L t +\ 4- t/* + i cosh (C49) 

2c,f+i 


Inserting equation (C48) into equation (C49) gives as a result 


Pt -\- 1 

U% — — sinh <p ( +iL ( +i 
ZcJ+l 

plz£i±i cqs ^ 

Z^ + 1 sinh <p ( +iL(+ { 


(' 


) 


cosh <p e+l L !+l 


(C50) 


In view of equation (C41), just upstream (left) of the cavity, 
the pressure and volume flow are given by 


1 


0 


. / ( M l \ 

\cjp a c\ k) 


Pi 


U e 


(C51) 


where U g is expressed by equation (C50). 

When the multiplication in equation (C51) is carried out, 
the pressure-velocity vector at the output of tube t is 


fa\yVt/ i\ 

j - — ( a + - I 

\c/p a C\ k) 


Pi 


IP, + ^ ** ^ 


^C,i+ 1 


Zcj+l ^f+A+1 


(C52) 


So we finally calculate the pressure at the input of tube i by 
multiplying equation (C52) by the ABCD matrix for a tube 
of length L t . The result is 


Pe~\ — p g cosh <p g L g + Z c ^ sinh <p g L g 

Pi + 1 


( /«\yF f / A p l+1 

X ) y ( ~ ) ( a + 7 )Pt + ^ 

/ \ c / Pa c \ kj Z c ,i + 1 


sinh <p g + \L ( + 1 


+ ( Pt Pl+l COSh <Pt+lLi+ ' 1 cosh <Pt+ ,Lf+ ! 


“) 


Z c f+1 sinh ipt+iLg+i 
Dividing by p t and rearranging yield 


(C53) 


( a + 1 )Z C> * sinh <p t L t 


Pi— i 

= cosh 4- / 

Pi V/PaC 


Z c ,i sinh f L r Pb-i \ 

~ — — [ cosh ) (C54) 

Z c ^ + i smh \ Pi 


(Cl 9) imply (when the numerator and denominator are 
multiplied by tube length L) 


<f___ n <pL p a c 1 

Z r — _ . 


■a 


r OJ 

y( 7 


y SL 


(C55) 


(As shown by the fact that it cancels out in equation (C55), 
the characteristic impedance Z c does not properly depend on 
the length of a tube. However, Bergh and Tijdeman apparently 
wish to introduce it so that they can express their results in 
terms of the tube volume V t = SL. Note that in all Bergh and 
Tijdeman expressions, the tube volumes V t are somewhere 
divided by the corresponding tube length L.) Finally, when 
equation (CIO) and the next-to-last part of equation (C19) are 
used, the ratio of characteristic impedances in the last term 
is equivalent to 


Zc,e ZiPi+x _ K,i+i<Pi+iL{F( a i+\) 

Z c ,t + 1 PfZ f+1 Vt,i<PiLt+ 


(C56) 


The second of the three terms summed on the right side of 
the equation can be transformed into the Bergh and Tijdeman 
form by noticing that the last equalities in equations (Cl 2) and 


So, regarding the F function as the ratio of Bessel functions 
we g et the Bergh and Tijdeman recursion relation of 
reference 4, which is 
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Appendix D 
Computer Programs 


Several short FORTRAN routines have been written to apply 
the ABCD matrices to the analysis of probe tubes coupled to 
transducer cavities. The basic package consists of three 
general-purpose subroutines, which are called by a main 
program designed to reflect the structure of the system being 
analyzed. Listings for the subroutines and two sample main 
programs are given at the end of this appendix. 

Subroutine TUBE evaluates the ABCD matrix of equation 
(C29) for a duct of generalized radius G, cross-sectional area 
S, and length L. It uses the formulas from either the circular- 
cylinder or the infinite-plate geometry, depending on whether 
the value of parameter IS is 1 or 2, respectively. The previous 
values of the components of the pressure-volume velocity 
vector are read in as the array PU, they are multiplied by the 
ABCD matrix, and the new vector is returned in the same 
variables. Information about the frequency, thermodynamic 
properties of the fluid, and mathematical constants is provided 
in COMMON blocks which must be initialized in the calling 
program. The propagation constant PHI calculated within 
TUBE is made available to the calling program in another 
COMMON block, to use if desired. 

Function RFUN is used to evaluate either the F of equation 
(B24) or the ratio J 2 fJo> depending on the value of parameter 
IS. It is called by subroutine TUBE and uses subroutines in 
the NASA Lewis Research Center FORTRAN library to 
compute the Bessel functions. For large arguments the 
common limiting form given in equation (21) is used. 

Subroutine VOL evaluates the ABCD matrix of equation 
(C41) for a small cavity in the line. The cavity is characterized 
by its volume V, ddd factor SIG, and poly tropic constant 
POLYK. The appropriate values of POLYK vary from 1 .0 
for an isothermal process to 1.4 for an adiabatic one. This 
constant is set equal to 1.0 in the programs of reference 5; 
perhaps a better estimate could be made on the basis of the 
expression for n in equation (B31). Subroutine VOL also takes 
in and updates the current pressure-flow vector. 

Two examples of calling programs are provided. The first 
was written to calculate, for frequencies from 25 to 5000 Hz, 
the pressure transfer function of a single tube and closed 
volume system, using both circular- and flat-geometry 
formulas for the tubing. The results are plotted on the same 
axes for comparison. The constants are defined in lines 1300 
through 2300, 2800 through 2900, and 3600 through 3800. 
Dimensions of the tubing are read in by lines 2950 through 
3190, and the transmission-line calculations are carried out 
at each frequency in the DO-loop in lines 3500 through 4960. 


The pressure-volume velocity vector for the circular-tubing 
model is the array PUR, while that for the flat-plate model 
is PUF. The remainder of the program plots the results by 
using routines standard at the Lewis Research Center. 
Numerical values for circular geometry computed with this 
program agree with calculations from the programs of 
reference 5 and with graphs in Bergh and Tij deman (ref. 4). 

The output of this program for sample sections of flat-oval 
tubing leading to a small, closed cavity is shown in figures 
9 through 12. In addition, to help an interested user determine 
that the computer listings have been copied correctly, table I 
gives selected numerical data calculated by this program. 

The second program is an example of how the routines 
would be used to compute the pressure response of an infinite- 
line probe. Once again the calculations are done twice, once 
for each model geometry, and both results are plotted for 
comparison. The sample graphs in figure 15 show how similar 
the results are. 


TABLE I. — SAMPLE NUMERICAL DATA FROM 
COMPUTER PROGRAM 


(a) Input Parameters 


Example 

Cavity 

volume, 

cm 3 

Generalized 
radius of 
tube, 
cm 

Area of 
tube, 
cm 2 

Length 
of tube, 
cm 

Frequency, 

Hz 

1 

0.01325 

0.060 

0.01714 

5.15 

1350 

2 

.01325 





5.15 

2700 

3 

.01325 





5.15 

4200 

4 

.01630 





7.71 

900 

5 

.01630 





7.71 

1900 

6 

.01630 


— 




7.71 

i 

2800 


(b) Calculated Output 


Example 

Circular tube 

Flat plates 

Amplitude 

ratio 

Relative 

phase, 

deg 

Amplitude 

ratio 

Relative 

phase, 

deg 

1 

8.98 

-91.1 

8.51 

-91.1 

2 

.905 

-179.2 

.904 

-179.2 

3 

4.46 

-261.7 

4.33 

-261.7 

4 

7.29 

-82.6 

6.82 

-82.6 

5 

.917 

-180.3 

.915 

-180.3 

6 

3.62 

-249.2 

3.49 

-249.2 
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(a) Ratio of amplitude of transducer cavity pressure to amplitude of inlet 
pressure. 

(b) Relative phase of transducer cavity pressure and inlet pressure. 

Figure 15.— Calculated ratio of cavity to inlet pressure and relative phase of two 
pressures as function of frequency for infinite-line probe made from flat-oval 
tubing with g = 0.061 cm and S = 0.0181 1 cm 2 . Distance from inlet end 
to transducer, 5 cm; distance from transducer to closed end, 100 cm. 


In general, names of variables are intended to suggest their 
significance. Definitions of variables placed in COMMON are 
as follows: 


AMDENS 

AMPRES 

FREQ 

GAMMA 

J 

J32 

OMEGA 

PHI 

RHOC 

SPEED 

SRPRN 

vise 

WAVENO 


ambient density, p a 

ambient pressure, p a 

ordinary frequency, / 

ratio of specific heats, y = c p /c v 

j raised to 3/2 power 
angular frequency, co 
propagation constant, <p 
the product p a c 
adiabatic speed of sound, c 
square root of Prandtl number, VPr 
coefficient of shear viscosity, p 
free-space wave number, k = co/c 
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Program Listings 


Subroutine 

0000100 
0000200 
0000300 
0 0 0 0 4 C 0 
0000500 
0000600 
0000700 
0000800 
0000900 
0001000 
0001100 
0001200 
0001400 
0001500 
0001600 
0001700 
0001800 
0001900 
0002000 
0002100 
0002200 
0002300 
0 0 0 24 Q 0 
0002500 
0002600 
0002700 
0002800 
0002900 
0003000 
0003100 
0003200 
0003300 
0 C 0 34 0 0 
0003450 
0003500 
0003600 
0003700 
0003800 
0003900 
0004000 
0004100 
0004200 
0004300 
0004400 
0004500 
0004600 
0004700 
0004800 
0004900 
0005000 
0005100 


TUBE — 


SUBROUTINE TUBEC G, S, L , IS , PU) 

C 

C SUBROUTINE TO EVALUATE ABCD MATRIX FOR TUBE OF 

C LENGTH L, CROSS SECTIONAL AREA S, AND 

C 2 ( AREA ) /PERIMETER EQUAL TO G. THE PRESSURE- 

C VOLUME FLOW VECTOR PU IS UPDATED. 

C 

C IS = 1 — > USE CYLINDRICAL DUCT FORMULAS 

C IS = 2 — > USE PARALLEL PLATE FORMULAS 

C 

COMPLEX PU(2),J,J32,PHI 
COMMON /XFREQ/ FREQ, OMEGA, UAVENO 

COMMON /THERMO/ GAMMA , SRPRN, AMPRES, AMDENS, SPEED, RHOC,VISC 
COMMON /XTUBE/ PHI 
COMMON /XCONST/ J,J32,PI 
C 

REAL L 

COMPLEX ABCD(2,2),PUTEMP(2),Z,ZC,AL,ALP,CN,CD,CE,CS 
COMPLEX RFUN 
C 

AL = J32X(G*SQRT (AMDENS* OMEGA/ VI SC) ) 

ALP = ALKSRPRN 

CN = GAMMA + (GAMMA - 1. )*RFUN(ALP, IS) 

CD = RFUN ( AL > IS ) 

PHI = WAVENOKCSQRTCCN/CD) 

C 

Z = -JXUAVEN0XRH0C/(SXCD) 

C 

C CHARACTERISTIC IMPEDANCE = SERIES IMPEDANCE/LENGTH / 

C PROPAGATION CONSTANT 

C 

ZC = Z/PHI 
C 

CE = CEXP(PHIXL) 

ABCDC 1,1) = (CE + CMPLX(1.,0.)/CE)/CMPLXC2.,0.) 

ABCD( 2,2) = ABCDC 1,1) 

CS = (CE - CMPLX(1.,Q.)/CE)/CMPLX(2.,0.) 

ABCDC 1,2) = ZCXCS 
ABCDC 2 , 1 ) = CS/ZC 
C 

DO 1 K = 1, 2 

1 PUTEMPCK) = PUCK) 

C 

DO 3 M = 1, 2 
PUCM) = CMPLXC 0 . > 0 , ) 

DO 2 N = 1, 2 

2 PUCM) = PUCM) + ABCDCM,N) X PUTEMPCN) 

3 CONTINUE 
C 

RETURN 

END 


Function RFUN — 


0 0 0 1 9 D 0 
CQ02000 C 
CQ02100 
0002200 C 
0002220 C 
0002240 C 
0002260 C 
0002265 
0 00227 0 
0002275 
0002280 
0002285 C 
0 002300 
0002400 C 
0 0 0250 0 
0 0 026 0 0 
0002700 C 
0002800 C 
0002900 C 
0003000 C 
0003100 
0003200 
0 0 0 330 0 
0 0 03400 
0 0 03500 
0 0 036 0 0 
0 0 0 37 0 0 
0 0 0390 0 
00 04 0 0 0 
0004100 C 
0004200 
0 0 0430 0 
0004500 
0004600 C 
0004700 
0 0 048 0 0 
0004900 C 
0 0 050 0 0 
0005100 C 
00 05200 
00 05300 


COMPLEX FUNCTION RFUNCZ,I5) 

COMPLEX Z,J0, Jl, J2,Y0,Y1,N,D,B 

IS = 1 “> EVALUATES J2(Z)/J0(Z) 

IS = 2 — > EVALUATES ( TAN (Z/2 )/ ( Z/2 ) ) - 1. 

A = REAL CZ) 

C = AIMAG(Z) 

IF (A .EQ. 0.0 .AND. C .EQ. 0.0) GO TO 3 
IF CADS(C) .GE. 140.) GO TO 4 

GO TO Cl, 2), IS 

1 CALL ZBE5JYCZ, JO, J1,Y0,Y1) 

J2 = CMPLXC2. 0,0. 0) x Jl/Z - JO 


SCALE NUMBERS BEFORE DIVIDING TO 
AVOID EXPONENT OVERFLOW 

ARO = REAL(JO) 

AIO = AIMAGCJO) 

AR2 = REAL ( J2 ) 

AI2 = AIMAGCJ2) 

A = AMAX1CABS(AR0),ABS(AI0),ABS(AR2),ABS(AI2)) 
N = CMPLXC AR2/A, AI2/A) 

D = CMPLXC ARO/A, AIO/A) 

RFUN = N/D 
GO TO 5 

2 B = Z/CMPLXC2. 0,0.0) 

RFUN = C5IN(B)/CB*CC05CB)) - CMPLXC 1 . 0 , 0 . 0 ) 

GO TO 5 

3 RFUN = CMPLXC 0.,0.) 

GO TO 5 

4 RFUN = CMPLXC '-1 . , 0 . ) + CMPLXC 0 2 . )/Z 

5 RETURN 
END 
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Subroutine VOL — 


00 0 0 1 00 
0000200 C 
0000300 C 
0000400 c 
0000500 C 
0000600 C 
0000700 C 
0000720 C 
0000740 C 
0000800 C 
0000900 
0001000 
0001200 
0001400 
0001500 C 
0001600 
0001700 C 
0001800 
0001900 
0002000 
0002100 
0002200 C 
0002300 
0002400 
0002500 C 
0002600 
0002700 
0002800 
0002900 
0003000 
0003100 C 
0003200 
0003300 


SUBROUTINE VOL ( V , SIG, POLYK , PU ) 

SUBROUTINE TO EVALUATE ABCD MATRIX FOR A LUMPED 
VOLUME V. PARAMETER SIG IS THE DIMENSIONLESS 
INCREASE IN TRANSDUCER VOLUME DUE TO DIAPHRAGM 
DEFLECTION WITH PRESSURE, AND PGLYK IS THE 
POLYTROPIC CONSTANT FOR THE VOLUME (SEE B & T) . 

THE PRES SURE- VOLUME VELOCITY VECTOR PU IS 
UPDATED , 

COMPLEX PUC2),J,J32 

COMMON /XFREQ/ FREQ , OMEGA , WAVENO 

COMMON /THERMO/ GAMMA , SRPRN , AMPRES , AMDENS , SPEED, RHOC, VISC 
COMMON /XC0N5T/ J,J32,PI 


COMPLEX ABCD(2,2),PUTEMPC2) 


ABCDC 1 , 1 ) 
ABCD ( 2 , 2 ) 
ABCDC 1,2) 
ABCDC 2,1) 


CMPLXC 1 . * 0 • ) 

ABCDC 1,1) 

CMPLXC 0 . , 0 . ) 

J*CWAVENO#GAMMA*V*CSIG + 1 . /PGLYK)/RHOC) 


DO 1 M = 1, 2 
1 PUTEMP CM ) = PU CM) 


DO 3 M = 1, 2 
PUCM) = CMPLXC 0 . , 0 . ) 

DO 2 N = 1, 2 

2 PUCM) = PUCM) + ABCDCM,N) x PUTEMPCN) 

3 CONTINUE 


RETURN 

END 


Program to plot comparison of circular and flattened tubes.— 


0000100 c 

0000200 C 
0000300 C 
0000450 C 
0000500 C 
0000600 
0000700 
0000900 
0001000 
0001100 
0001200 C 
0001300 
0 0 0 1400 
0001500 
0001600 C 
0001700 
0001800 
0001900 
0001925 
0001950 
0 0 020 0 0 
0002100 
00 02200 
0002300 
0002400 C 
0002500 
0002520 
0002530 
0002540 
0002560 
0002580 
0002600 C 
0002300 
0002900 
0002930 C 
0002950 
0002970 
0002990 
0003010 
0003030 
0003050 
0003070 
0003090 
0003110 
0003130 
0003150 
0003170 
0003190 
0003310 C 
0003420 
0003440 
0003460 C 
0003500 
0 00 360 0 
0003650 
0003700 
0 G 0 38 0 0 
0003900 C 
0004000 
0004100 
0004200 
0004220 
0004240 
0004300 
0004350 
0004400 C 
0004500 
0004600 
0004650 
0004700 


PROGRAM TO PLOT COMPARISON OF CIRCULAR AND 
FLATTENED TUBES , CALCULATIONS USE 
TRANSMISSION LINE FORMALISM. 

CGS UNITS USED THROUGHOUT. 

COMPLEX J,J32,PHI 

COMMON /XFREQ/ FREQ , OMEGA , WAVENO 

COMMON /THERMO/ GAMMA , SRPRN , AMPRES, AMDENS, SPEED, RHOC, VISC 
COMMON /XTUBE/ PHI 
COMMON /XCONST/ J,J32,PI 

PI = 3 . 141593 
J = CMPLXC 0 . , 1 , ) 

J32 = CEXPC CMPLXC 0. ,0.75*PI)) 

GAMMA = 1.4017 
SRPRN = 0.8414 
ATM = 1.0129E6 
AMPRES = ATM 
TEMPC = 26.85 

AMDENS = AMPRES/C2 , 8688E6XCTEMPC+273 . 15) ) 

SPEED = SQRT ( GAMMA* AMPRES/AMDENS ) 

RHOC = AMDENS * SPEED 
VISC = 1.846E-4 

COMPLEX PRATIO, PURC2),PUFC2) 

REAL XFC200),ARC20Q),AFC200),PR(200),PF(200) 

INTEGER I VARS (10) 

REAL XTITC3)/’ FREQ’ , ’UENC’ , ’Y ’/ 

REAL YTIT1(4)/’AMPL’,’ITUD’,’E RA’,'TIO ’/ 

REAL YTIT2C2)/'PHAS’ , 1 E V 

SIG = 0. 

POLYK = 1.2 

3 WRITE (6, 1) 

1 FORMAT C/ f ENTER CAVITY VOLUME 1 ) 

READ (5, 2, END=20) V 

2 FORMAT CG10) 

WRITE (6, 4) 

4 FORMAT (/' ENTER GENERALIZED RADIUS OF DUCT 1 ) 

READ (5, 2, END=2Q) G 

WRITE (6, 6) 

6 FORMAT C/’ ENTER CROSS-SECTIONAL AREA OF DUCT’) 

READ (5, 2, END=20) S 

WRITE (6, 7) 

7 FORMAT C/’ ENTER LENGTH OF DUCT 1 ) 

READ (5, 2, END=20) XLEN 

RMX = 0. 

NPTS = 200 

DO 10 KF = 1, NPTS 
FREQ = FLOAT ( 25*KF) 

XF(KF) = FREQ 

OMEGA = 2. * PI * FREQ 

WAVENO = OMEGA/SPEED 

PURC1) = CMPLXC 1 . , 0 . ) 

PURC2) = CMPLXC 0 . , 0 . ) 

CALL VOL CV , 5IG, POLYK, PUR) 

PUF(l) = PURC1) 

PUFC2) = PURC2) 

CALL TUBE(G,S, XLEN, 1, PUR) 

CALL TUBE(G,S,XLEN,2,PUF) 

PRATIO = CMPLXCl. ,0. J/PURCl) 

ARCKF) = CABS ( PRATIO ) 

RMX = AMAX1CRMX, ARCKF) ) 

PHA = ATAN2(AIMAG(PRATIQ),REAL(PRATIO))X180./PI 
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0 0 0480 0 

0004850 

0004900 C 

0 0 0 4 9 1 C 

0004920 

0004930 

0004940 

0004950 

0004960 

0004970 C 

0004973 

0004976 

0004979 

0004982 

0004985 

0004988 

0004991 

0004994 

0004997 

0005000 

0005003 

0005006 C 

0005009 

0005012 

0005015 

0005018 C 

0005021 

0005024 

0005027 

0005030 

0005033 

0005036 C 

0005039 

0005042 

0005045 

0005048 

0005051 

0005054 C 

0005057 

0005060 

0005063 

0005066 C 

0005069 

0005072 

0005075 

0005200 C 

0005220 

0005240 C 

0005260 

0005300 

0005400 


IF (PHA . GT . 0.) PHA = PHA - 360. 

PR(KF) = PHA 

PRATIO = CMPLXCl. >0.)/PUF(l) 

AF(KF) = CABS ( PRATIO ) 

RMX = AMAXKRMXi AF(KF) ) 

PHA = ATAN2 (AIMAG( PRATIO ), REAL (PRATIO) ) *180 ./PI 
IF (PHA .GT. 0.) PHA = PHA - 360. 

10 PF(KF) = PHA 

I VARS ( I ) = 8 
I VARS ( 2 ) = NPTS 
I VARS ( 3 ) = 66 
IVAR5C4) = 62 
I VARS ( 5 ) = NPTS - 1 
I VARS ( 6 ) = 25 
I VARS (7 ) = 0 
IVARSC8) = NPTS 

CALL GINTVL(0. , RMX, 10,1, AMIN, AMAX) 

CALL GINTVL(0., 5000., 10,0, AMIN, AMAX) 

CALL GPLOT(XF,AR,IVARS) 

I VARS ( 3 ) = 98 
I VARS ( 4 ) = 65 
CALL GPLOT(XF,AF,IVARS) 

CALL TITLE(4,9,15,XTIT) 

CALL TITLE(3,15,15,YTIT1) 

CALL CGRNER(l) 

CALL COPY < 1 ) 

CALL DISPLA(l) 

I VARS ( 3 ) = 66 
I VARS ( 4 ) = 62 

CALL GINTVL(-360. ,0. ,10,1, AMIN, AMAX) 

CALL GINTVLCO. ,5000. >10,0, AMIN, AMAX) 

CALL GPLOTCXF, PR, I VARS) 

I VARS ( 3 ) = 98 
I VARS ( 4 ) = 65 
CALL GPLOT (XF, PF, I VARS) 

CALL TITLE(4,9,15,XTIT) 

CALL TITLE(3,5,15,YTIT2) 

CALL DISPLA(l) 

GO TO 3 

20 CALL TERM 
STOP 
END 


Program to plot response of infinite line systems — 


0000100 c 
0 0 0 0 ?. 0 0 c 

0000300 C 
0 0 0 0400 
0 0 0 050 0 
0 0 0 06 0 0 
000 0700 
00 0 0800 
0000900 C 
0 0 0 1 000 
0 0 0 1 1 0 0 
0 0 0 120 0 
0001300 C 
0 0 01400 
0 0 0 150 0 
0 0 0 1 6 0 0 
0 0 0 17 00 
0001800 
0 0 0 190 0 
0 0 020 0 0 
0002100 
0002200 
0002300 C 
0002400 
0002450 
0002500 
0002600 
0002700 
0002800 
0002900 
0003000 C 
0003100 
0003200 
0003300 C 
0003400 
0003500 
0003600 
0003700 
0003800 
0003900 
0004000 
0004100 
0004200 
0004300 
0004400 
0004500 
0004600 
0004700 C 
0004800 
0004900 
0005000 C 
0005100 
0005200 
0005300 
0005400 
0005500 
0005600 C 
0005700 
0005800 
0005900 
0006000 
0006100 C 
0006200 
0006300 
0006400 C 
0006500 
0006600 
0006700 C 
0006800 
0006900 


PROGRAM TO PLOT RESPONSE OF INFINITE LINE 
SYSTEMS. 

COMPLEX J,J32,PHI 

COMMON /XFREQ/ FREQ , OMEGA , WAVENO 

COMMON /THERMO/ GAMMA , SRPRN , AMPRES , AMDENS , SPEED, RHOC , VISC 
COMMON /XTUBE/ PHI 
COMMON /XCQNST/ J,J32,PI 

PI = 3.141593 
J □ CMP LX ( 0 . ,1.) 

J32 = CEXP(CMPLX(0.,0.75#PI)) 

GAMMA = 1.4017 
SRPRN = 0.8414 
ATM = 1.0129E6 
AMPRES = ATM 
TEMPO = 26.85 

AMDENS = AMPRES/ ( 2 . 8688E6* (TEMPC+273 . 15 ) ) 

SPEED = SQRT ( GAMMA X AMP RES/ AMDENS ) 

RHOC = AMDENS * SPEED 
VISC = 1.846E-4 

COMPLEX PRATIO,PURC2) ,PUF(2) 

COMPLEX PREFR,PREFF 

REAL XFC20 0 ) , AR ( 20 0 ) , AFC 20 0 ) , PRC 20 0 ) , PFC 20 0 ) 

INTEGER I VARS (10) 

REAL XTIT(3)/ f FREQ f , 'UENC' ,'Y V 

REAL YTIT1(4)/ , AMPL F aiTUD 1 , T E RAS'TIG V 

REAL YTIT2C2)/ f PHAS' ,'E V 

SIG = 0. 

POLYK = 1.2 

3 WRITE (6, 1) 

1 FORMAT C/ f ENTER DISTANCE FROM SENSOR TO OUTLET 1 ) 

READ C5, 2, END=2Q) XL1 

2 FORMAT (G10) 

WRITE (6, 4) 

4 FORMAT (/' ENTER GENERALIZED RADIUS OF DUCT 1 ) 

READ (5, 2, END=20) G 

WRITE (6, 6) 

6 FORMAT (/' ENTER CROSS-SECTIONAL AREA OF DUCT 1 ) 

READ (5, 2, END=20) S 

WRITE (6, 7) 

7 FORMAT C/ f ENTER LENGTH OF INFINITE TUBE 1 ) 

READ C5, 2, END = 20) XL 2 

RMX = 0. 

NPTS = 200 

DO 10 KF = 1, NPTS 
FREQ = FL0ATC25XKF) 

XF (KF) = FREQ 

OMEGA = 2. * PI * FREQ 

WAVENO = OMEGA/SPEED 

PUR ( 1 ) = CMPLXC1.0E-6, 0. ) 

PURC2) = CMPLX(0.,0.) 

PUFC 1 ) = PUR ( 1 ) 

PUFC2) = PUR ( 2) 

CALL TUBE(G,S,XL2»1,PUR) 

CALL TUBE(G,S,XL2,2,PUF) 

PREFR = PUR(l) 

PREFF = PUFC 1 ) 

CALL TUBE(G,S,XL1,1,PUR) 

CALL TUBE(G,S,XL1,2,PUF) 
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0007000 C 
0007100 
0007200 
0007300 
0007400 
0007500 
00076 00 
0007700 C 
0007800 
0007900 
0008000 
0008100 
0008200 
0008300 
0008400 C 
0008500 
0008600 
0008700 
0008800 
0008900 
0009000 
0009100 
0009200 
0009300 
0009400 
0009500 
0009600 C 
0009700 
0009800 
0009900 
0010000 C 
0010100 
0010200 
0010300 
0010400 
0010500 
0010600 C 
0010700 
0010800 
0010900 
0011000 
0011100 
0011200 C 
0011300 
0011400 
0011500 
0011600 C 
0011700 
0011800 
0011900 
0012000 C 
0012100 
0012200 C 
0012300 
0012400 
0012500 


PRATIO = PREFR/PURC1) 

ARCKF) = CABS ( PRATIO ) 

RMX = AMAX1CRMX, ARCKF) ) 

PHA = ATAN2(AIMAG(PRATIO),REAL(PRATIO))#180./PI 
IF (PHA .GT. 0.) PHA = PHA - 360. 

PRC KF) = PHA 

PRATIO = PREFF/PUFC 1 ) 

AF(KF) = CABS (PRATIO) 

RMX = AMAXHRMX, AF(KF) ) 

PHA = AT AN2 ( A IMAGC PRATIO) * REAL ( PRATIO ) )#180 ./PI 
IF (PHA .GT. 0.) PHA = PHA - 360. 

10 PF(KF) = PHA 


IVARS 

I VARS 

IVARS 

IVARS 

IVARS 

IVARS 

IVARS 

IVARS 

CALL 

CALL 

CALL 


( 1 ) = 8 

(2) = NPTS 

(3) = 66 

(4) = 62 

(5) = NPTS - 1 

(6) = 25 

(7) = 0 

(8) = NPTS 
GINTVLCO., RMX, 10,1 
GINTVLCO, >5000. ,10 
GPLOTCXF, AR, IVARS) 


AMIN, AMAX) 

0, AMIN, AMAX) 


IVARSC3) = 98 
IVARSC4) = 65 
CALL GPLOT(XF,AF, IVARS) 

CALL TITLE(4,9,15,XTIT) 
CALL TITLE(3,15,15,YTIT1) 
CALL CORNER(l) 

CALL COPY(l) 

CALL DISPLAC1) 


IVARSC3) = 66 
IVARSC4) = 62 

CALL GIN TV L (-360. ,0. ,10, 1, AMIN, AM AX) 
CALL GINTVLCO. ,5000. ,10,0, AMIN, AMAX) 
CALL GPLOTCXF, PR, IVARS) 

IVARSC3) = 98 
IVARS ( 4 ) = 65 
CALL GPLOTCXF, PF, IVARS) 

CALL TITLE(4,9,15,XTIT) 

CALL TITLEC3,5,15,YTIT2) 

CALL DISPLA(l) 


GO TO 3 


20 CALL TERM 
STOP 
END 
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